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ABSTRACT 

New improvements to compute Mie scattering quantities are presented. They are 
based on a detailed analysis of the various sources of error in Mie computations and on 
mathematical justifications. The algorithm developed on these improvements proves to 
be reliable and efficient, without size (x = 2nR/X) nor refractive index (m = mR — imi) 
limitations, and the user has a choice to fix in advance the desired precision in the results. 
It also includes a new and efficient method to initiate the downward recurrences of Bessel 
functions. 
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1. INTRODUCTION 



The Mie theory of light scattering by a homogeneous sphere is used for many prob- 
lems of atmospheric optics and also in other fields in Physics. The application of Mie 
theory still needs modern computers for numerical calculations of the many functions and 
coefficients involved. The primary difficulty is in the precise evaluation of expansion coef- 
ficients a n and b n . This is further aggravated as x gets large, and when the calculation of 
size distribution is needed. An optimization of computer time for reliable computation is 
clearly of necessity. 

The formulas for Mie scattering are well known 1 ' 2 . Here we follow the notation of 
Bohren and Huffman 3 . The scattering and extinction efficiency factors are given by 

2 N 

^ = ^E( 2n + 1 Hl a -l 2 + l 6 -l 2 ) 

x n=1 (i) 

2 N 

Qe = — J^C 271 + X ) Re ^ + ^ 

X 

n=l 

where x = 2nR/\ is the size parameter of the problem, R being the radius of the sphere, A 
the wavelength of the light and N a large enough number. The Mie scattering coefficients 
a n and b n are functions of x and the relative refractive index m = tur — wii, with 
rriR > 1, ™>i > 0. 

x( n (x)ip' n (y) - y(' n {x)^ n {y) 
b = y^n(x)^' n (y) - xip' n (x)il) n (y) 

y(n(x)lp'n(y) - xC^{x)lp n {y) 

where y = mx and ^ n (z) , £ n (z) are the Riccati-Bessel functions related to the spherical 
Bessel functions j n (z) and y n (z): 

(n(z) = zj n (z) - izy n (z) 

These functions are known in closed form (Ref. 4, p. 437) but it is more convenient to use 
the recurrence relation 

X n+1 (z) = F n {z)X n {z) - X n _!(z), 
F n {z) = (2n + l)/z . 

where X is any of the functions in eqn. (3). 

Presently, there are many versions of Mie scattering computer codes (Dave 5 ' 6 , 
Blattner 7 , Grehan and Gouesbet 8 ' 9 , Wiscombe 10 ' 11 , Goedecke et al. 12 , Miller 13 ) and au- 
thors who had been doing Mie calculations (Kattawar and Plass 14 , Deirmendjian 15 , Quen- 
zel and Muller 16 , Bohren and Huffman 3 ). These are reflected in performing our work. 
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Some essential points should be addressed by any Mie scattering algorithm: 

1) How to determine the number iV for truncating a Mie series. 

2) Whether the Riccati-Bessel functions will be computed by upward recursion or by 
downward recursion. 

3) If downward recursion is used, how to initialize it. 

4) How to structure the algorithm in an efficient way. 

Answers to all the above questions constitute the objective of this paper. We focus 
particularly on analyzing the numerical error sources and show that our Mie algorithm 
permits users to prescribe a precision e beforehand, to effect an efficient, reliable Mie 
coefficients calculation. Needless to say, the precisely evaluated Mie coefficients a n , b n are 
required for calculating the angular scattering amplitudes 1 ' 2 ' 3 ' 5 ' 6 ' 10 . 

2. CONVERGENCE PROPERTIES OF THE MIE SERIES 

In this section we shall estimate the error introduced in some typical quantity such 
as the efficiency factors, by keeping a finite number iV of partial waves in the Mie series. 
We shall also find a criterion for choosing the value of N. In this section the quantities 
a n , b n themselves are assumed to be computed exactly. 

In order to investigate the convergence properties of the scattering coefficients a n , b n 
we shall make use of very well known properties of the spherical Bessel functions (e.g. ref. 
4, p. 438 and ff.). Let us recall some properties which are relevant for us: 

i) 

lim ip n {z) =0, lim ( n (z) = oo . (5) 

ii) For z = x real, ip n ( x ) an d Cn( x ) have two distinct regimes as functions of n: 

a) oscillating regime for n < x. ip n ( x ) and ( n ( x ) keep changing their sign regularly, 
and |^ n (x)| and iCnO^O! are bounded by slowly changing functions of n. 

b) exponential regime for n > x. ip n ( x ) becomes exponentially decreasing and ICnO^)! 
becomes exponentially increasing. 

In view of these considerations one concludes from eqn. (2), that all the partial 
waves n < x (x being the size parameter from now on) will contribute to the Mie series 
and convergence will appear only after n enters in the exponential regime. This is so 
because ip n (x), ip' n (x) go very quickly to zero in the numerator and (n( x ),C n ( x ) g° to 
infinity in the denominator. On the other hand ipn(y) > ip'niy) a PP ear both in numerator 
and denominator and therefore seem to play no role in the convergence. We can emphasize 
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this fact by writing 



_ j) n (x) r , V'n(^) - s/y) + - yA n (x) 

Cn(x) n ( n (x) n(y/x-x/y) + xA n (y)-yB n (x) 

__ ipn{x) , , _ ip n {x) yA n (y) - xA n (x) 
Cn(x) n ( n (x) yA n (y) - xB n (x) 

where we have extracted the factor ip n (x) / Cn( x ) responsible for the convergence of a n and 
b n and also we have reexpressed the ratios tp' n (z) /ip n (z) and C n ( x )/Cn(x) in terms of (ref. 
4, p. 439) 



Let us state more clearly our assumption: we shall assume that the quantities 
[a] n , [b] n are bounded by slowly varying functions of n in the exponential regime n > x. 
The validity of this assumption will be analyzed in a later section. 

If [a] n and [b] n are well behaved for large n, we can approximate them by their 
asymptotic values in order to discuss the convergence of a n and b n . In order to take ad- 
vantage of this approximation we can use the asymptotic expansion of the Bessel functions 
for large orders (ref. 4, p. 365), 

A n (z) ~ F n (z), B n (x) ~ F-\x) as n ^ oo (8) 

where the next term in the expansion has a higher power of 1/n. We obtain 

[a] n ~i^ + 0(-), [6] B ~0(-) (9) 

In practice, for x < n < N , [a] n and [6] n are both of the order of unity, (unless m is nearly 

1, in which case [a] n , [b] n ~ 0). On the other hand, recalling that wir > 1, it can be proved 
1 1 — 2 1 

that | - l _ l _ m 2 1 < 2, therefore a good enough estimate is 



[a] n , [&]„«! (10) 



Using this and the asymptotic values (8), it can be shown that the truncation error in Q e 
is bounded by 

SQ e < \a N \ (11) 

The proof is presented in Appendix I where it is shown that the series YI^Ln+i l a ™l 
converges faster than some geometric series. Let us note that what actually appears in Q e 
is Rea n , not \a n \, therefore the bound (11) will usually be conservative. This is especially 
true for small m; because in this case Rea n ~ |a n | 2 (i.e. Q e ~ Q s ) and |a n | 2 \a n \ for 
n > N. 
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Let us now find a criterion for choosing the number N of partial waves that should 
be taken into account. For this purpose let e be the error allowed in the calculation, and 
let us take 5Q e as a typical quantity in the problem. Then N should be taken so that 



SQ e < € (12) 

Taking the quantity Q e has the advantage of being simple and also that 5Q S < 5Q e , 
because |a n | 2 < Rea n (i.e. Q s < Q e for each partial wave). Other interesting quantities, 
such as the scattering amplitudes, have similar convergence properties as Q e and Q s . 

Putting together the bound (11), the criterion (12) and the estimate (10) we find 
the following prescription 



Cn(x) 



< c (13) 



In order to find something more convenient let us make use of the Wronskian identity (ref. 
4, p.439) 

tpn(x)(n-l(x) - 1pn-l(x)Cn(x) = i, (14) 

and the asymptotic values of A n (x) and B n (x). In this way we obtain (within approxima- 
tions keeping the order of magnitude) 

MxKnW^-iF-^x), (15) 

This allows us to remove ip n (x) from (13) and finally we obtain the prescription for N 

llmOvOr)] > \[\, (16) 

which has been written in a form convenient for being checked while ( n (x) is being com- 
puted by upward recurrence. In getting (16) we have neglected a factor Fn(x) from (15) be- 
cause by doing so iV may increase at most by one unit (recall that ( n (x) j Xn-i(x) ~ F n (x)). 
Also we have used that Ke( n (x) = ifj n (x) is negligible as compared to lm( n (x) in the ex- 
ponential region. 

It is remarkable that the value of iV obtained from (16) for e = 10 -8 is virtually 
identical to the standard prescription N = x + ex 1 / 3 + 1, with c = 4.3. It is shown in 
Appendix II that it must be so using asymptotic expansions for Cn{x), and also how to 
modify c if some other precision e is desired. To know N(x) in advance is necessary if the 
computer code is to be vectorized 10 ' 11 . 

3. NUMERICAL ERROR AND UPWARD RECURRENCE 

In this section we shall discuss the propagation of numerical error through the 
calculation. 
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It is known that the determination of ip n (z) by upward recursion is intrinsically 
unstable (see e.g. ref. 5). Let us clarify this point.* For the sake of simplicity let us 
assume that the numerical error is coming from the initial values 

ipo(z) = i/) (z) + e , 4>i{z) = ip!{z) + ei (17) 

but the recursion itself is free of roundoff error, i.e. 

ip n+ i(z) = F n (z)4> n (z) - ^n-l(z) (18) 

eo, 6i being small numbers depending on the precision of the computer, and ip n (z) being the 
numerical sequence that is actually obtained instead of the exact one, ip n ( z )- Subtracting 
the exact recursion for ip n (z) from (18) we find 

5ijj n+1 (z) = F n (z)Sip n (z) - (fyn-i (*) (19) 

where 8ip n (z) = tp n (z) — ip n (z) is the error in our numerical sequence. Any sequence 
satisfying the recurrence relation (4) is a linear combination of ip n (z) and ( n (z), therefore 

5iP n (z)=r 1 Mz)+v'Uz) (20) 

The small numbers rj, rj' are directly related to eo, e\ through eqn. (17), namely 

V = KeoCi - dCo) ^ 21 ^ 
rj' = -i(e V>i - eiV'o) 

Recalling now that ( n (z) diverges for large n we conclude that the absolute error in ifj n (z) 
will eventually blow up. More generally, if the recursion itself is not exact due to computer 
roundoff error, ifj n (z) is rather given by 

^n{z) = VnMz) + VnCn(z) (22) 

where rj n , rj' n are of the order of the roundoff error or the initial values error, whichever the 
largest. In any case the conclusion is still that Sip n (z) is small for small n (or while n is in 
the oscillating regime for z nearly real), but blows up when n enters in the exponentially 
increasing regime of (n(z). Since ifj n (z) itself goes to zero in the exponential regime, tp n (z) 
has less and less correct figures at each step. 

We can extract some corollaries from the previous discussion: 

1) The upward recursion is always unstable for computing ip n (z) for large n, de- 
pending on z. The error 8ip n (z) grows as |£ n (z)|. On the other hand the upward recursion 
is perfectly stable for computing ( n (x) for any value of n. This is because S( n (x) still 

* We thank one of the referees for providing us with a simpler proof of this statement. 
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grows as |£ n (a;)|, therefore the relative error in ( n (x) is kept small. Note however that the 
relative error in the quantity ReCn(^) = i>n( x ) is n °t at all small. 

2) A downward recursion is stable for computing ip n (z), because |Cn(-z)| is either 
slowly changing (in the oscillating regime) or quickly decreasing with decreasing n (in the 
exponential regime) . This allows for taking even very rough estimates for the initial values 
of ip n (z) in the downward recursion and the ratio ip n -i(z)/ipn ( z ) w ih still quickly approach 
the exact value A n (z). On the other hand, a downward recursion is not appropriate for 
computing ( n (x) or the ratio B n (x) if it starts in the exponential regime. 

Now let us study the influence of the numerical error on the a n , b n coefficients, and 
hence on Q e if an upward recursion is used to compute ip n (x). In this analysis ( n (x) 
and B n {x) are assumed to be exact due to previous considerations. On the other hand 
A n {y) is also assumed to be exact. The effect of using approximate values of ip n (y) will be 
considered later. We can make the discussion for a n . Similar conclusions will hold for b n . 
Eqn. (6) can be rewritten as 

a n = ^-f(A n (x)), (23) 

where only the A n (x) dependence is shown explicitly as it is the only relevant one for error 
analysis. The relative error in a n will be given by 

~ ^IL + iL^AlL (24) 

a n Ipn f A n 

Recalling the definition (7), the relative error in A n can be estimated to be of the same 
order of magnitude as that of ip n , and taking into account that / is a smooth function of 
the order of unity (cf. eqn. (10)), one gets the estimate 

Sa n « a n —- « a n r] — = rj f w rj . (25) 

where use has been made of eqn. (22) and rj' is some typical value of rj' n . 

This means that the absolute error in a n or b n , remains roughly constant throughout 
the computation. Of course eqn. (24) holds only for small 5ip n: but this is guaranteed 
as iV is of the order of x and so the recurrence does not go deep inside the exponential 
region. The important consequence of eqn. (25) is that the upward recursion can be used 
to obtain tp n (x) because the error introduced is of the order of the roundoff error (see 
however the comment at the end of Section 6). Let us note that this fact is consistent 
with available algorithms for doing Mie calculations, where ip n (x) and ( n (x) are always 
computed by upward recursion (e.g. refs. 5,11). 

Let us consider now the effect of the numerical error coming form tp n (y). We have 
argued before that an upward recursion would not be appropriate for computing i/j n (z) in 
general, however we have just shown that it can be used in the case of ip n (x). The reason 
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for this was that the relative error in ip n (x) grew as C n (x)/tfj n (x) but the quantities a n and 
b n themselves converged to zero as ip n ( x ) / Cn( x ) • Both factors cancel rendering 8a n and 
5b n bounded. We cannot apply a similar argument to 5ip n (y) and therefore an upward 
recursion is not reliable to compute ip n (y) for arbitrary y. We can consider two limiting 

cases 

a) mi = 0. In this case y is real and greater than x, thus the instability in ^ n {y) starts 
only after that in ifj n (x), therefore the upward recursion can be used. 

b) Large m/. From the initial values 4 

ipo(z) = sin(z) , tpi{z) = - sin(z) — cos(z) 



z 



Co{z) = iexp(-iz) , d(z) = ( ^ - 1 ) exp(-iz) 



(26) 



one can see that ip n ~ exp(mix), ( n ~ exp(— mix), for small n, thus ifj n is much 
larger than ( n . On the other hand eo,i are related to the computer precision, 
typically eo,i ~ h/>o,i with r pa 10~ 16 in double precision. Upon substitution in (21) 
we find that r\ is small but r\' ~ rexp(2m/a;) which is not necessarily small. The 
relative error in ifj n (z) goes as 



5lp n {z) 
1pn(z) 



^o(z) 



Co(z) 



Cn(z) 



1pn(z) 



(27) 



For small n the relative error is small, of the order of r, however for n ~ \z\, where 
ifj n and ( n are of the order of unity, the relative error is r^o/Col ~ rexp(2m/x) 
which is large for large mi. Therefore the upward recursion is not stable in this 
case. 

To summarize, the upward recursion to compute ifj n (y) can be used if mi is small 
enough but becomes unstable for large mi. We have not analyzed in any detail in which 
cases the upward recursion for ifj n (y) is reliable, therefore we shall only consider downward 
recurrences for this quantity. See however refs. 10,11 for an extensive analysis of this 
problem through computer experiments. Noting that all we need is the ratio A n (y), for 
1 < n < N, we can use the downward recursion 

A n (y) = F n (y) - — J-p- . (28) 

Computing the initial value A]^{y) requires some algorithm such as that of Lentz 17 or the 
one we present in the next section. Let us estimate now the precision required in A^{y) in 
order not to introduce an error in Q e larger than the prescribed precision e. By arguments 
similar to those used for ip n (x), we have 

5a n 5A n (y) , . 

— w . , / 29 
a n A n (y) 
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where 5a n is the error introduced by SA n (y). Given that the downward recursion is stable 
we can assume that 



5a 7 



< 



SA N (y) 



A N (y) 



for n<N 



(30) 



Using this relationship one gets for the numerical error in Q e 



N 



SQ e ~ — 5^(2n + l)Sa n < Q t 



n=l 



SA N (y) 



A N (y) 



(31) 



Therefore the numerical error from A n (y) will be under control by imposing 



SA N (y) 



A N (y) 



< 



(32) 



Let us note that this criterion will be conservative in general. An exception would be the 
case of y being real and bigger than N. In this case the recurrence (28) has no healing 
properties (for it already starts in the oscillatory regime) and hence the equal sign is 
reached in (30). 

4. INITIALIZATION OF THE DOWNWARD RECURRENCE 

In this section we present a new method to compute An(z), of similar efficiency to 
that due to Lentz 17 (actually ours needs one multiplication less at each step). This method 
has the advantage of being able to implement a precision condition as that in eqn. (32), 
hence controlling the required precision in A n (y). 

Let X n (z) and Y n (z) be two sequences satisfying the recurrence (4) for some value 
of z (the dependence on z is irrelevant here) . Then they will satisfy the Wronskian identity 



C — X n Y n+ i — X n+ iY n 



(33) 



where C is independent of n. We can rewrite it as a difference equation 



C — YnY-r 



n 1 n+1 



n+1 



(34) 



and solve it in X r 



X n = DY n + CY n J2(YkY k+1 )- 1 : 



(35) 



D being a constant. To write (35) we have assumed that Y n is a sequence going to infinity 
for large n, which is true for almost any solution of the recurrence (4). If we take Y n as a 
fixed sequence and regard C, D as free parameters, then X n is the most general solution of 
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the recurrence relation (4). In particular for D = 0, X n goes to zero as n goes to infinity, 
as a consequence it must be proportional to tp n , 



Mz) = C(z)Y n (z) ^(y fc (z)y fc+1 (^))- 1 (36) 

fc=n 

The constant C cancels after computing the ratio A n (z) 

oo 

A n (z)=Y- 1 {Y n . 1 +Y- 1 [J2(y k Y k+1 )- 1 }- 1 }. (37) 

k=n 

Finally, a simpler formula can be obtained for An(z) by choosing as starting values for 
the sequence Y n 

Y N _ X = , Y N = 1 (38) 

oo 

A N (z) = [J2(yk(z)Y k+1 (z))- 1 ]- 1 - (39) 

k=N 

About the convergence of the series in (39), we note that it is very fast when Y k enters 
in its exponential regime. Note that for real y the convergence begins only after k > y. 
A similar conclusion was reached by other authors 11 in Lentz's method which basically 
follows the same principle as ours and so has similar convergence properties. 

The sequence in eqn. (39) must be truncated at some value k = M in such a 
way as to fulfill the requirement (32). This can be easily done by noting that the error 
introduced in A^ 1 (y) is of the order of the last term taken into account (this follows from 
\Y k /Y k -i\ ~ \F k \ >2 for large /c), 

SArf » (y M ^M+i) _1 . (40) 
On the other hand we should require 

-!SA N 



\SA- N \y)\ 



A 



N A 



N 



< 



F N (y) Q t 



(41) 



where we have made use of eqns. (8) and (32). Recall now that for x > 1, F N and Q e are 
of the order of unity whereas for x <^ 1 the product of F^Qe is still of the order of unity, 
therefore the final criterion to truncate (39) is 

\(Y M (y)Y M+1 (y))- 1 \ < e. (42) 

To finish this section we shall show how to avoid ill-conditioning in (39), which will appear 
if Yfc gets too near to zero for some value of k. To do this we can use the recurrence relation 
(4) to write 

i + _!_ = n-i + n + i = f„ 



Yk-iY k Y k Y k+ i Y k -iY k Y k+1 Y k -iY k+ i ' 
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which is well behaved even for Y k = 0. 

5. COMPUTATIONAL ALGORITHM 

Using the previous ideas, we have developed a computational algorithm which we 
shall briefly describe now. The input is x, m and e and the main output are the coefficients 
a n and b n , and N. To start with, analytic expressions for Co( x ) an d CiO^) are taken 
to initiate an upward recurrence for ( n (x)- This quantity is kept in a (complex) array 
variable. The recurrence stops when the condition (16) is fulfilled, providing the value 
of N. The quantities ip n (x) are automatically obtained as the real part of ( n (x)- As a 
second step, A N (y) is computed using eqns. (38), (39) and (42). Here we note that from a 
computational point of view an equivalent form of (42) is more convenient, which consist 
in doing the check for the absolute values of the real and imaginary parts. This is much 
faster than computing the modulus of a complex number. 

Then a downward recurrence is performed for A n (y), eqn. (28), until n = 1. Si- 
multaneously, a n and b n are computed using ( n (x) and A n (y). The quantities Q s and Q e 
can then be computed. We have not developed any especial algorithm for computing the 
scattering amplitudes Si and S 2 . To do this efficiently see for instance ref. 11. 

The criteria developed above are intended to be robust, hence they are rather con- 
servative. As a consequence the error in Q e is smaller than the prescribed precision e. This 
is especially true for small values of x, whereas for x ^> 1, about two more figures than 
expected are obtained. We point out also that Q s is always obtained as accurately as Q e 
or more. This fact was expected because the criteria were stated for \a n \ while Q s goes as 
|a n | 2 which converges faster. 

6. RESONANT TERMS IN THE MIE SERIES 

Let us recall that after eqn. (7) we stated a smoothness assumption for the quantities 
[a] n , [b] n , namely that they are nearly constant in the x exponential regime and do not 
play any role in the convergence of the Mie series, which is only controlled by the ratio 
ip n ( x ) / Cn( x ) • In particular this assumption implied that the highest partial wave with a 
relevant contribution is independent of m (cf. eqn. (16)). In other words, N is a function of 
x only. This result is also supported numerically, (see for instance refs. 10,11). Therefore 
it was a surprise for us to discover that strictly speaking such a statement must be false. 
Moreover, for any choice of N as a function of x only, and for any prescribed value of n, 
n > N ', one can always pick a value of m (in fact infinitely many of them) in such a way 
that the n-th term in the Mie series is not negligible, for instance one can make a n = 1. 
The consequence of this that in order to guarantee that the numerical value of Q e is correct 
within some prescribed precision, N should depend on m as well as on x. 

In order to clarify the point let us consider the worst case, which is also the simplest, 
namely m; = 0, i.e. y real. This is the only case in which \a n \ or \b n \ can reach the value 
1. The point can be made for a n : recalling that for z real Re( n (z) = ip n (z), eqn. (2) can 
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be rewritten as 



a r 



Re D T 

D„. 



D n = x( n (x)ip' n (y) - yC n (x)ip n (y) 
where D n is a complex quantity. Obviously a n = 1 if and only if 



(44) 



lmD r 







(45) 



Let us regard x and n as given and look for solutions of (45) in the variable y. The equation 
can be rewritten as 

l<(y) _ HmCQr) (46) 

y^n(y) xlm( n (x) 

In the interval y > n, ip n (y) is a rea *l oscillating function of y with infinitely many zeroes. 
Between two zeroes of ip n (y), the l.h.s. of eqn. (46) takes every real value, therefore there 
are infinitely many solutions to our equation for any values of x and n, no matter how large 
is n as compared to x. For these values of x, m, and n, a n will not at all be negligible. 

Let us now show that these resonances do not occur for unrealistic values of m. 
Typically (and asymptotically for large y) the distance between two consecutive zeroes of 
ip n (y) is of the order of n, therefore for given x and n the lowest resonant value of m will 
occur near the interval (^, n ^ 2L ) approximately. For large x this happens for m near to 
unity, and all the other resonant values will follow at a distance of about ir/x from each 
other. 

From a rigorous point of view these findings would invalidate the estimates (10) and 
their consequences. They would also invalidate any algorithm in which N depends on x 
only, namely every existent algorithm known to us. In fact the only practical way to make 
sure that the resonant partial waves have been accounted for would be to take N greater 
than y in order to guarantee that ^ n {y) has no zeroes for n > N. 

Nevertheless it is clear that in practice the existent algorithms to do Mie scattering 
calculations work. To account for this fact we should consider not only the existence of 
resonant partial waves but also their width. Let us show that for sensible choices of N (as 
a function of x) and for n > N the resonances are so narrow that they will not normally 
show up. Let y be one the values of y such that a n = 1. A look to eqn. (44) shows 
that for generic y, ReD n goes as ^ n {x) whereas D n goes as ( n (x), therefore a n is very 
small. However for the especial value yo there is a cancellation between two huge numbers 
in ImD n , leaving a n of the order of unity. The range of values of y for which a partial 
cancellation takes place is related to the slope of D n in y = yo, namely 



(47) 



D n 




ReD n 




tpn(x) 


D'n 


y=yo 


D' n 


y=yo 


Cn{x) 



Where D' n = dD n /dy. In other words, if N is large enough only by a very careful choice of 
m or x can one find one these resonant contributions. More precisely, recalling eqn. (13), 
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we can see that m or x should be fine tuned at least with a precision e in order to pick 
a resonant term for some n > N . On the other hand, except for these rare cases, a n , b n 
are indeed small and of the order of ip n (x) / Xni x ) ? therefore our analysis applies. If m is 
allowed to be complex, a more involved analysis would be needed, but we expect that the 
conclusion would not differ. 

Let us finally note another consequence of the resonant terms on the calculation, 
even when they are taken into account. For one of these terms the quantity / in eqn. (23) 
is no longer of the order of unity, on the contrary it is rather large, and the last step in 
eqn. (25) cannot be taken. This means that a resonant term amplifies the error due to 
ifj n (x). The cure is simply to compute ^ n (x) by downward recursion for x < n < N. This 
has in fact been observed in selected quantities such as the backscattering efficiency for 
suitable values of x and m (Ref. 5). 

7. CONCLUSIONS 

In this paper we have addressed several points relevant to Mie scattering calcula- 
tions. To be specific: 

a) We have estimated the error introduced in the calculation by truncating the Mie 
series, thereby finding a prescription for choosing N. We have found that in the generic 
case N depends on x only. 

b) The possible instabilities in the recursions used to compute ip n and £ n have been 
analyzed. We have found that upward recursion is always unstable for computing ifj n (z) 
if n is large enough. However it can be used to compute ipn(x) in Mie calculations. As a 
matter of fact ip n (x) is computed in this way in nowadays available algorithms. We have 
also found that upward recursion can be used for ifj n (y) if m; is small enough, but no 
criterion is given for how small m; should be. 

c) A criterion has been established for the allowed error in ifj n -i(y)/ifj n (y). 

d) A new method to compute tfj n -i(y)/^ n (y) is presented which is efficient and 
allows for controlling the error and removing ill-conditioning. 

e) It has been shown the existence of resonant terms in the Mie series which can also 
appear for n > N. Strictly speaking the existence of these terms invalidates any algorithm 
in which A is a function of x only. However we have also shown that those resonant terms 
are extremely rare, namely they appear with a probability of the order of e. 

A specific algorithm is also described. It is meant to be robust and efficient for a 
wide range of size parameters and refractive indices. With this algorithm we have written 
the computer program LVEC-MIE 18 , which is available both in single and double precision 
contacting V.E. Cachorro. 
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APPENDIX I 



Let us justify the bound (11). To do so we shall study the convergence rate of the 
terms left out in the series, n> N. In this region we can make use of the estimate (10), 



$Qe = — Yl (2^ + l)Re(a n + 6 n ) 

X n=N+l 

2 °° 

< ^ E (2n + l)(|a n | + |6„|) 



n=N+l 

oo 



2>2 ^ > 

n=N+l 



II 



8 



X' 



E 

n=N+l 



n 



Cn(x) 

B n (x) B n -i{x) B N+1 (x) 



A n (x) A n _ 1 (x) ' ' ' A N+1 (x) 



Cn(x) 



(1.1) 



Now making use of (8) and recalling that F n (x) is a monotonically increasing function of 
n, we obtain 



E 



1 1 



n=N+ 
x. 

n=iV+l 
"I — A 

X 



n 



^At(x) 



Civ(x) 



2(N-n) 



^jv(ic) 



+ 



2 \F 2 N (x)-l (F*(x)-l) 



(1.2) 



For small x, N = 1 and F N (x) is large, hence 



5(5 e < 2|a r 



(1.3) 



on the other hand, for large x, N xs x and F N (x) ~ 2, 



5Q e < ^-|a n | 
o a; 



(1.4) 



In both cases eqn. (11) is valid (up to factors of the order of unity). 
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APPENDIX II 



In order to know in advance the value of iV that will be obtained from the prescrip- 
tion (16) for given x and e, let us recall that Im( n (x) = a/ nx/2 Y N+ i(x), Y v {z) being the 
Bessel function of the second kind. Let v and c be defined by 

N = v- - 

2 (ILl) 
x = v - cz/ 1/3 . 



Note that for large v, eqn. (II. 1) can be inverted to give iV rs x + cx 1 ^ 3 . Now we can make 
use of the leading order term in the asymptotic expansion of Y v for large v and fixed c, 
ref. 4, p. 367: 

1/6 

(II.2) 



Im6v(aO~-v^(T) Bi(2 1 / 3 c) , 



where Bi(z) is the Airy function of the second kind, ref. 4, p. 446. This function is given 
by 

Bi(z) = z- 1 ^f(z)eMl^ 2 ) ' (IL3) 
where f(z) is nearly constant for z > 1 with f(z) xs l/y/n, ref. 4, p. 449. Thus 

1/6 



Im Cat (#) 



2 w 



(II.4) 



The right hand side of (II. 4) has a very strong dependence on c whereas it depends very 
smoothly on v. Actually (z//2y / 2) 1 / 6 is of the order of unity for v = 1 up to 10 5 . Therefore 
using eqn. (16), c will be determined by e. We find that c = 4.3 corresponds to e = 10 -8 . 
Other values are c = 4.0, e = 10~ 7 , and c = 5.0, e = 10~ 10 , computed for v = 100 in 
(II.4). 
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